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Abstract 

The richness of natural images makes the quest for optimal representations in image pro- 
cessing and computer vision challenging. The latter observation has not prevented the design 
x-/ of image representations, which trade off between efficiency and complexity, while achieving 

C/3 accurate rendering of smooth regions as well as reproducing faithful contours and textures. The 

, most recent ones, proposed in the past decade, share an hybrid heritage highlighting the multi- 

scale and oriented nature of edges and patterns in images. This paper presents a panorama of 
the aforementioned literature on decompositions in multiscale, multi-orientation bases or dictio- 
naries. They typically exhibit redundancy to improve sparsity in the transformed domain and 
sometimes its invariance with respect to simple geometric deformations (translation, rotation). 
Oriented multiscale dictionaries extend traditional wavelet processing and may offer rotation 
invariance. Highly redundant dictionaries require specific algorithms to simplify the search for 
an efficient (sparse) representation. We also discuss the extension of multiscale geometric de- 
compositions to non-Euclidean domains such as the sphere or arbitrary meshed surfaces. The 
etymology of panorama suggests an overview, based on a choice of partially overlapping "pic- 
tures". We hope that this paper will contribute to the appreciation and apprehension of a 
stream of current research directions in image understanding. 
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1 Introduction: Vision Aspects, Scope and Notations 

1.1 Background on Vision Aspects of Scale 

Many natural-world object features are substantive only over a certain spatial extent. In other 
words, the scale of observation is crucial in object recognition and understanding. For instance, a 
chair would be easily recognizable in the scale of a few meters. But neither at a centimeter scale 
which captures the chair's texture and not its object appearance, or at a hectometer scale, where 
the chair's appearance is hardly distinguished from other surrounding objects. 

Accordingly, early neurophysiological studies in biologic perception reveal that those objects 
are generally apprehended differently according to the scale of observation by the sensory receptors 
and the cortex of mammalians [HE]- Efficient information extraction is thus required for artificial 
sensing systems to mimic standard biologic tasks such as object recognition. 

Pixel-based representations as linear combinations of "delta" functions suffice for simple data 
manipulation but are very limited for higher level tasks. Only assuming some sufficient resolution 
in the data, the lack of prior knowledge in the extent of objects to be analyzed calls for tools 
able to unveil the appropriate scales and to allow a hierarchical representation of the underlying 
features [31 EJ E]. Disregarding the peculiar fractal formalism [6l [7] where similar phenomena 
appear at different scales (what is called self -similarity), special attention has been paid to data 
transformations able to capture object features over a range of scales in a more compact form. 
Sparsity, amounting to a reduced number of parameters in a suitable domain, is thus used as a 
heuristic guide to image understanding. Bearing analogies with findings in vision processes jS], 
several sparse decompositions have proven efficient in image compression, with the discrete wavelet 
transform (DWT) as their most well-known avatar, often intermingled with information theory and 
technical wizardry, from bit plane arithmetic coding [9] to trellis coded quantization. A compact 
history and a paper collection are given in [TQ1 [TTj , respectively. 

Yet, beyond image compression transforms, other decomposition techniques are needed, with 
more resolving power in complex scene detection, denoising, segmentation or, in a broad sense, 
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(a) 



(b) 



Figure 1: Two faces of the cartoon-texture model: (a) Yogi bear (b) Fingerprint. 

scene understanding. As a matter of fact, standard separable wavelet transforms appropriately 
detect point-like (0-D) singularities and address mild noise levels. Still they generally lack perfor- 
mance in dealing with higher dimensional features combining both regularity and singularity such 
as edges, contours or regular textures, that may also be anisotropic. Amongst their limitations are 
shift sensitivity, limited orientation selectivity, rigid and uneven atom shapes (e.g., fractal-looking 
asymmetric Daubechies wavelets), crude frequency direction selection. Major challenges reside in 
a proper definition of the underlying regularity (with respect to each feature) and corresponding 
singularities. These challenges are amplified by additional degradations from which acquired data 
may suffer such as blur, jitter and noise. Descriptive mathematical models of images combining 
cartoon and textures become increasingly popular jT2j 13J and progressively yield tractable algo- 
rithms. We note that there exists a continuum of real- world images between cartoon and textures, 
ranging from cartoon-ish Yogi bear in Fig. |l(a)| to "texturaP fingerprints in Fig. |l(b)[ In between 
these two extreme image types, there exists many possible variations in image object complexity. 

Moreover, both contours and (even regular) textures are known to be ill-defined. They are 
indeed viewer- and scale-dependent concepts in discrete images or volumes. Consider an image 
resulting from a combination of piecewise smooth components, contours, geometrical textures and 
noise. Their discrimination is required for high level image processing tasks. Each of these four 
components could be detected, described and modeled by different formalisms: smooth curves or 
polynomials, oriented regularized derivatives, discrete geometry, parametric curve detectors (such 
as the Hough transform), mathematical morphology, local frequency estimators, optical flow ap- 
proaches, smoothed random models, etc. They have progressively influenced the hybridization of 
standard multiscale transforms towards more geometric and sparser representations of such compo- 
nents, with improved localization, orientation sensitivity, frequency selectivity or noise robustness. 

1.2 Scope of the Paper 

Geometry driven "*-let" transforms [14] have been popular in the past decade, with a seminal 
ancestor in [15]. Early [16J, a debate opened on the relative strength of Eulerian (non-adaptive) 
versus Lagrangian (adaptive) representation, now pursued with the growing interest in dictionary 
learning [17J. 

As of today, the authors believe that the discussion is not fully settled in the various different 
uses of sparsity in images. Neither has the trade-off between redundancy and sparsity. A number 
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of early papers on geometric multiscale methods appear in [18J. Comparisons are drawn in fT9ll20]. 
while [2TJ E21 ESI E21 Ell f°cus on ridgelets, curvelets and wedgelets, as representative of fixed 
and adaptive decompositions. The present paper aims at providing a broader panorama of the 
recent developments in multiscale decompositions targeted to efficient representation of geometric 
features in images: smooth content (multiscale or hierarchical), edges and contours (locally spatial) 
and textures (locally spectral). We emphasize the main characteristics and differences pertaining 
to spatial, directional and frequency selectivity of the selected methods. The paper therefore 
cites a dense set of references, ranging from continuous to discrete representations, from (nearly) 
orthogonal to (fully) redundant. As a guiding thread to this panorama, we illustrate some of the 
reviewed geometric multiscale decompositions on a memorial plaque 1 in Szeged University, Hungary, 
depicted in Fig. [2j It features simple objects (embedded rectangles and a disk), a few differently 
oriented features and regular textures at different scales. Since some of the illustrations have been 
slightly enhanced to improve the clarity of details, they are available in original resolution online 
[26] . This picture finally honors Alfred Haar's originative paper [27J Zur Theory der orthogalen 
Funktionen Systeme (On the Theory of Orthogonal Function Systems) and his eponymous wavelet. 
He also founded Acta Scientiarum Mathematicarum together with Frigyes Riesz, whose works 
percolated wavelet theory 
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Figure 2: Szeged University Memorial plaque in honor of A. Haar and F. Riesz: A szegedi matem- 
atikai iskola vildghiru megalapitoi (The world-wide famous founders of the mathematical school in 
Szeged). 



The paper is organized as follows: the remaining of Section[T]is devoted to context and notations 
for image representations. Then, as a preliminary to geometric tools, a quick survey of early 
multiscale decompositions is presented in Section [2| More recent transforms, termed directional 
or geometrical, circumventing aforementioned drawbacks, are discussed in Section [SJ Owing to 



1 Courtesy of Professor Karoly Szatmary, http : / / astro . u- szeged . hu/ szatmary . html who performed scalograms 
analysis of variable stars as early as in 1992 [25] . 
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the additional degrees of freedom provided by these representations, a discussion is carried out in 
Section [4] on redundancy and adaptivity. The extension of frequency, scale and directionality to 
non-Euclidean spaces or grids such as the sphere, are presented in Section [5| Finally, concluding 
remarks are given in Section [6j 

1.3 Mathematical Framework 

1.3.1 Notations and Conventions 

This paper describes numerous mathematical methods designed for different spaces and ge- 
ometries. We have tried therefore to adopt coherent representations for the many mathematical 
notions that coexist in this text. For instance, functions and vectors in high dimensional spaces 
are generally referring to some signal of interest (e.g., 1-D signals or images). They must therefore 
share the same notations and we thus decided to write them as simple lowercase Roman or Greek 
letters. However, coordinate systems, vectors in 2 or 3 dimensions and multi-indices are denoted 
in bold symbols. 

The (Hilbert) space JL 2 (X) is the space of square integrable functions on the space X, i.e., given 
the (Lebesgue) integration measure dp on that space, L, 2 (X) = h 2 (X,dp) = {/ : X —> C : \\f\\ 2 := 
$x \f( u )\ 2 &P( U ) < °°}- I n L 2 (Af) the inner product between two functions g, h G ~h 2 (X) is denoted 
by (g, h) = j x g*(u) h(u) dp(u) with * the complex conjugation. By extension, for p ^ 1, we also 
use the (Banach) spaces 1?(X) = U>(X,dp) = {f : X -> C : := J x \f(u)\ p dp(u) < oc}, with 

We also use some discrete spaces as the common £ P N = (C N , || • with \\v\\p := J2i \ v i\ P f° r 
p ^ 1 and v G C N , with again the shorthand ||-|| — ||-||2. In £ 2 N , the inner product between u, v G l 2 N 
is written (u, v) = u • v = J2 u i v i- Whether the overused notations (•,•) or || • \\ p are applied to 
continuous or discrete mathematical objects will remain clear from the context. The spaces t v are 
the generalization of the previous finite spaces to infinite sequences, i.e., £ p = {v = (viji^N '• \\ v \\p — 

£»^oN < °°}- 

For functions / G h 2 (X) or discrete sequences v G l 2 N , f and v denote the Fourier transform 
of / or v respectively. For instance, for X = R and / G L 2 (R), f(u) = -j= J R /(*) e~ lujt dt 

and f(t) = -j= J R f(uj)e lujt duo are the Forward and Inverse Fourier transform respectively. For 
/ G L 2 (R 2 ) and x,uj G R 2 , the same transforms are f(u>) = ^ J R2 f(x)e~ luJ ' x d 2 x and f(x) = 
^ f$2 f(u)e lu ' x d 2 u?. For v G the same transforms are = ^7= Vj exp(— 2'Kijk/N) and 
Vj = ^= J2k exp(27ri jk/N). In matrix algebra notations, this can be rewritten as v = Tv and 
v = T*v, where the Fourier matrix T G C NxN is given by T^k — exp(27ri jk/N), and J 7 * is its 
complex adjoint. The convolution by time-invariant filter h operates as (f*h)(t) = f(u)h*(t — 
u)du and (v ★ h) n — h* n ,v n _ n i in continuous and discrete sample domain 2 respectively. The 
ubiquitous Gaussian kernel with scale parameter a > is denoted by G a (x) — exp(— ^2 ||x|| 2 ), 
with G(x) = Gi(x). 

2 With periodization for finite length vectors. 
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1.3.2 Image Representations in Bases and Frames 

Stability and Frames This paper describes processing methods that make use of a decomposi- 
tion of the image / E L 2 ([0, l] 2 ) into a family of atoms B — {^ m } m . Each atom ^ m E L 2 ([0, l] 2 ) is 
parameterized by a multi-index m (that might take into account its frequency position, scale and 
orientation). Numerical processing is performed on discretized images which are vectors / E M. N , 
where N stands for the number of pixels. The atoms of B are also discretized and the continuous 
inner products are replaced by the standard discrete inner product in M. N . 

To guarantee a stable reconstruction from the coefficients {(^m, /)}m 5 the family B is assumed 
to be a frame [29l EDI EHl ED |32j of L 2 ([0, l] 2 ) or R N , which means that there exist two constants 
< /xi ^ /i2 < oo such that for all / 

mll/ll 2 ^£l(Vv„,/>| 2 ^HI/ll 2 - (i) 

m 

Atoms are allowed to be linearly dependent, thus corresponding to a redundant representa- 
tion. Redundancy enables atoms to meet certain additional constraints, for instance smoothness, 
symmetry and invariance to translation or rotation. 



Thresholding for Approximation and Processing Using a dual frame {V>m}m [2H1, an image 
is recovered from the set of coefficients as / = /)^ m . The computation of the set of 

coefficients {(^ m , f)}m f° r a discrete image / E M. N is usually performed using a fast algorithm, 
that also enables a fast reconstruction of an image from coefficients. 

The basic processing operation, used in denoising and compression applications, is the thresh- 
olding 

f M = H T (f, B)= &™ /) $™ ( 2 ) 

m : \(ipmJ)\>T 

where M = # {m : |(^ m , /)| > T} counts the number of non-zero coefficients in Q. 

When ii\ = /^2, the frame is said to be tight (Parseval tight frame). If furthermore [i\ — fi2 — 1, 
then one can choose ^ m = and B — {^ m } m is then an orthonormal basis if \\ipm\\ = 1 for 
all m. In this last case, B performs the least energy reconstruction of Jm in ([2]), or equivalently, 
/m is the best M-terms approximation of /. The decay of the approximation error ||/ — /m\\ is 
related to both the average risk of a denoiser, and the distortion rate decay of a coder, see for 
instance [33J. This motivates the search for bases or frames B which can efficiently approximate 
large classes of (natural) images. When the frame is redundant, more complicated decomposition 
methods improve the sparsity of the representation (see Sec. 4.1). 



2 Early Scale-Related Representations 

2.1 Frequency, Heat Kernel and Scale-Space Formalism 

At the heart of modern signal processing techniques is the concept of signal representation, i.e., 
the selection of an efficient "point of view" in the study of signal properties that is not restricted 
to straightforward spatial descriptions. 

The most obvious alternative signal representation is its frequency reading, i.e., the one provided 



by the Fourier transform of the signal explained in Sec. 1.3.1 [3U|35]. However, this representation 
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is not sufficiently "local". It is indeed rather difficult to detect what spatial part of an image 
contributes to high peaks in the Fourier spectrum. Fig. [3] represents the amplitude spectrum 3 of 
the luminance component from Fig. [2| It exhibits a mixture of prominent vertical and horizontal 
directions with tiny fuzzy diagonal ones. 



Figure 3: Magnitude of the 2-D Fourier transform of the Haar-Riesz Memorial plaque in Fig. [2j 

An approach for obtaining a better localization is to introduce a notion of "scale" in the image 
observation. This has been performed very early in image and signal processing by either windowing 
or introducing scales in the Fourier transform |36jl3T] or observing a well-known diffusion process like 
the heat dynamics governed by the famous Heat equation. The idea relies on considering the image 
as an initial configuration of heat that is diffused with a time variable r > and in interpreting this 
time parameter as the "scale". Indeed, in this dynamic diffusion, small image structures will be 
smoothed early at small evolution time while larger ones persist for a larger duration. Interestingly, 
this diffusion is equivalently described by a filtering process: the convolution of the image by a 
Gaussian function G a of width a = y/2r [38l EH HQ]. This image unfolding into a scale-space 
domain has led to many new image processing techniques such as edge, ridge and feature detection 
[lU H2] . This is illustrated in Fig. [IJ where the original image is convolved with three different 
Gaussian kernels in dyadic progression. Large objects such as the white rectangular plaques persist 
across all scales, while brick and grid textures vanish in Fig. 4(c)| The overall redundancy of the 



Gaussian pyramid is given by the number of smoothing kernels. Taking advantage of the resolution 
loss, the redundancy factor may be reduced by sub-sampling, leading to the "Gaussian pyramid" 
construction. 

The scale content of the image can be decomposed further by computing, for instance, differences 
between two filterings performed at two different scales. This led to the famous Littlewood-Paley de- 
composition, or to the (invertible) Laplacian pyramid conveniently combining multiple sub-sampled 
low-pass filterings of images, creating a pyramidal scale hierarchy [43J. Interestingly, the resulting 
decomposition represented in Fig. [5] is a complete image representation that can advantageously be 
processed before reconstructing a new "restored" image (e.g., in image denoising). Additionally, 
image singularities are enhanced at fine scales, with low activity regions associated with coefficients 

3 The original image has been multiplied by a 2-D raised-cosine type apodizing window in order to reduce border 
discontinuity effects. 
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(a) (b) (c) 

Figure 4: Gaussian scale-space decomposition of the Haar-Riesz Memorial plaque at three different 
scales. 

being close to zero. Fast implementations of deformable (steerable or scalable) decompositions [H] 
are available for instance with recursive filters [H] or efficient multirate filter banks Wf \ H5 1 H9] . 




Figure 5: Laplacian pyramid decomposition of the Haar-Riesz Memorial plaque. 

Remarkably, the notion of Scale-Space has been defined and "axiomatized" more than 50 years 
ago by the Japanese mathematicians Iijima and Otsu, as presented in [50J. As we will realize 
throughout this paper, this scale-space representation (refer to [51J for a recent overview and 
axiomatic generalization) was the starting point of many new ways to represent images. 



2.2 Isotropic Continuous Wavelet Transform 

The continuous wavelet transform somehow generalizes the previous scale-space formalism 
driven by the Gaussian kernel to any "function" with enough regularity. The continuous wavelet 
transform was initially developed for the transformation of 1-D signals [52J and further extended in 
2-D first with isotropic wavelets. The case of non-isotropic (directional) wavelets was defined later 
[53] ( see Sec.lEO). 
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(a) (b) 

Figure 6: (a) The Marr wavelet (or Mexican hat), (b) Marr Wavelet singularity detector of the 
Haar-Riesz Memorial plaque. 

In one dimension, a wavelet ip is an integrable and well-localized function of L 2 (R), generally 
described as locally oscillating, i.e., J R ip(t)dt = 0. It may be dilated or contracted by a scale factor 
a > and translated to a position b E R: V>(&,a)(^) = 7^ V>(^r)- 

The continuous wavelet transform of a signal / E L 2 (R) probes its content with a "lens" ip^^) 
of zoom factor a and location b. Mathematically, 

W f (b,a) = [ f(t) j= ri^dt = (V( M ,/). (3) 

Interestingly, provided that -0 is admissible, i.e., when the two constants = 27T J^°° ^^^-dcj < 
oc are finite and equal 4 , that is, — — < oc, the signal / may be recovered from the 
coefficients Wf(b, a): 

r+oo r 

= i / / ^HM ( 4 ) 

This integral representation involves wavelets at every location and all positive dilations, i.e., / is 
decomposed on the continuous set of functions {ip(b,a) a £ Rfj b E R}. Many different kinds of 
(admissible) wavelets may be selected. We may cite the derivatives of Gaussian (DoG), the Morlet 
and the Cauchy wavelets, etc. Their selection is driven by the features to be elucidated in the data, 



e.g., frequency content with the Morlet wavelet or singularities with DoGs (Fig. 6(a) ) as illustrated 5 
in Fig. |gfb)l 

In two dimensions, the most natural extension of the 1-D-CWT is obtained by considering 
isotropic wavelets, i.e., wavelets r/j E L 2 (R 2 ) such that ip(x) = ^ r ad(||^||)> with x = (^1,^2), 
for some radial function ^ ra d : R + — » R. In that case, the wavelet family is generated by 2-D 
dilations and translations, i.e., we work with i/j(b,a)( x ) — a^(^T^) ^ na ^ are co Pi es °f V 7 translated 
to b = (h,b 2 ) E R 2 and dilated by a > 0. The 2-D CWT of the image / E L 2 (R 2 ) is then simply 



4 When ip is sufficiently regular, this condition reduces to a zero-average requirement, that is, f R ip(t) dt 
5 The YAWTb toolbox has been used, see http://rhea.tele.ucl.ac.be/yawtb/ 
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Wf(b, a) = (ip(b,a)i f) an d reconstruction of / is guaranteed by 

p+oc p 

/(*) = $/ / W>(6,a)^ (6ia) (x)d 2 6^, (5) 



if = (27t) 2 / r2 |^(^)| 2 /||^|| 2 d 2 fc < oc. The isotropic CWT is a useful analysis tool for edge 

2 I 



detection in images. For instance, by taking the (admissible) Marr Wavelet ip(x) — A[exp — i||#" 2 



(with A the 2-D Laplacian) also called Laplacian of Gaussian or Mexican Hat (see Fig. |6(a) ), the 



CWT of an image / acts as a multiscale edge detector. The topic of 1-D and 2-D continuous 
wavelet transforms is covered in more details in [52l |53j [5H l55| [33]. 

2.3 Discrete Scale-Space Representations 

Numerical computation requires that continuous expansions such as ^ and ^ be discretized. 
In this section, we detail some parameter samplings, such as dyadic or translation invariant grids. 
Together with a suitable choice of the wavelet function, they lead to stable representations where 
the original signal can be perfectly reconstructed from its coefficients. 

2.3.1 Multiresolution Analysis (MRA) 

In the context of a dyadic sampling where a — 2 3 and b = n2 J for j, nEZ, the canonical way 
to design a suitable wavelet function ij) in 1-D makes use of a multi-resolution analysis (MRA). 
It is defined as a nested sequence of closed vector subspaces (Vj)j eZ in L 2 (R) verifying standard 
properties [56]. Multiresolution analysis of a signal / consists of successively projecting the signal 
onto subspaces Vj in a series of increasingly coarser approximations as j grows. The difference 
between two successive approximations represents detail information. It amounts to the information 
loss between two consecutive scales, which lies in the subspace Wj, the orthogonal complement of 
Vj in Vj-i such that: 

Vj-^VjQWj. 

Then, with additional stability properties, there exists a wavelet ij) E L 2 (R) such that B — 
{2 _J / 2 '0(2 _J x — n) : n E N} is an orthonormal basis for Wj. 

2.3.2 Separable Orthogonal Wavelets 

A 2-D orthogonal wavelet basis B = {ipm}m of L 2 (R 2 ) for m = (j, n, k) is parameterized by a 
scale 6 2 J (j E Z), a translation 2 J n = 2 J (724,712) (n E Z 2 ) and one of three possible orientations 
k E {V, i7, D} 5 loosely denoting the vertical, horizontal and (bi) diagonal directions, the latter being 
poorly representative. Wavelet atoms are defined by dyadic scalings and translations ^ m (cc) = 
x — n) of three tensor-product 2-D wavelets 

^ V (x) = i/j(xi)</)(x2), ^ H {x) = (j){xi)^{x 2 ), and ^ D {x) = ^(xi)^(x 2 ), 

where (j) and are respectively 1-D orthogonal scaling and wavelet functions, see [5U EH [33] . When 
the scale interval is limited to j < J for some J E Z, the basis B is completed by the functional 
set A = {0(j, n )}n, with the 2-D separable scaling function <j)(x) = (j>{x{)(j)(x2). This set gathers all 



Here and throughout the rest of the paper, we use the convention that scale increases with j, as in s = 2 J . The 
converse convention is also often used in the literature. 
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the coarse scale wavelet atoms with j ^ J. The standard cascade image is depicted in Fig. [7| It is 
now critically sampled, i.e., free from redundancy (compare Fig. |5| and |6(b)| ). The approximation 
coefficients in A, a coarse image approximation at scale J, are represented in the bottom-left square 
of Fig. [7[ The other squares in this picture, associated to the "bands" {V, H, D} for j < J, exhibit 
some sparsity (few important coefficients), and horizontal and vertical edges are relatively well 
captured. 




Figure 7: Dyadic wavelet decomposition of the Haar-Riesz Memorial plaque. 

A non-linear approximation /m = Ht(/, B) in an orthogonal separable wavelet basis is efficient 
for smooth images or images with point- wise singularities. The approximation of a piecewise 
smooth image with edges of finite length decays like ||/ — /m|| 2 = 0(M~ 1 ). This result extends to 
functions with bounded variations [58] . and is asymptotically optimal. This decay is nevertheless 
not improved when the edges are smooth curves, because of the fixed ratio between the horizontal 
and the vertical sizes of the orthogonal wavelet support. 

2.3.3 Fast Algorithms for Finite Images 

A finite discretized image / e C NlxN2 of N = N ± N 2 pixels fits into the MRA framework 
by assuming that the pixel values of f n on n = (711,712) are the coefficients (0(j }7l ),/) of some 
continuous function / E L 2 (R 2 ) at a fixed resolution Vj, where 2 2J = N. 

The coefficients (^ nJ /) of / for j ^ J are computed from the discrete image / alone. This com- 
putation is performed using a cascade of filters interleaved with downsampling operators [56J . For 
compactly supported wavelets, this requires O(N) operations. Symmetric bi-orthogonal wavelet 
bases with compact support ease the implementation of non-periodic boundary conditions [59] . 
For infinite impulse response (IIR) wavelet niters, computations in the Fourier domain require 
0(N\og(N)) operations [60J, while recursive implementations |61| allow signal-adaptive implemen- 
tation. 
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While separable wavelets are not optimal for approximating generic edges, they lie at the heart 
of early state-of-the-art methods for compression and denoising. The JPEG 2000 coding standard 
[62] performs an embedded quantization of wavelet coefficients, and uses an adaptive entropic 
coding scheme that takes into account the local dependencies across wavelet coefficients. The 
sub-optimality of wavelets for the sparse representation of edges can be alleviated using block 
thresholding of groups of wavelet coefficients [63], that gives improvements over scalar thresholding. 
Advanced statistical modeling of wavelet coefficients leads to denoising methods close to the state- 
of-the-art, see for instance [6 ^ 165 } 166]. 



2.3.4 Translation Invariant Wavelets 

Given a discrete frame B — {^ m } m of C^, B is translation invariant if — r) E B for any 
-0 E B and any integer translation r. This property tends to reduce artifacts in image restoration 
problems like denoising, since, for such invariant frame, the thresholding operator Ht(/, B) becomes 
itself translation invariant. Discrete orthogonal wavelet bases described in the previous sections are 
not translation invariant and many authors have worked on recovering this useful capability. 

For instance, cycle spinning, proposed by Coifman and Donoho in [67], reduces wavelet artifacts 
by averaging the denoising result of all possible translates of the image, thus resulting in a transla- 
tion invariant processing. For an orthogonal basis B = {^ m } m , this is equivalent to considering a 
tight frame which is the union of all translated bases {'0 m (' — ^)}m,r- For a generic basis, this frame 
has up to TV 2 atoms. For a wavelet basis, the frame has 0(N\og(N)) atoms, and the coefficients 
are computed with the fast "a trous" algorithm in 0(Nlog(N)) [601 EH] . The translation invariant 
paradigm additionally draws a connection between the scale-space formalism (Sec. |2.1[) [69] and 



thresholding (Sec. 1.3.2). Several 2-D design described in the next sections attempt to (approxi- 



mately) address invariance (translation/rotation) without sacrificing computational efficiency. 



3 Oriented and Geometrical Multiscale Representations 

The variety of oriented and geometric multiscale representations proposed over the last few 
years requires broad grouping, arranged as follows: Sec. |3.1| presents directional methods closely 
related to 1-D decompositions. In Sec. |3.2| the directionality is addressed with diverse non-separable 
schemes. Finally, in Sec. |3.3[ directionality is attained by an anisotropic scaling of the atoms that 
yields various efficient edge and curve representations. 



3.1 Directional Outcrops from Separable Representations 
3.1.1 Improved Separable Selectivity by Relaxing Constraints 

As discussed in Sec. |2.3.1| discrete orthogonal wavelets may be viewed as a peculiar instance 
of orthogonal filter banks [70]. A well-known limitation in 1-D is that orthogonality (hence non- 
redundant), realness, symmetry and finite support properties cannot coexist with pairs of low- and 
high-pass filters, except for the Haar wavelet. 

We decide to briefly mention here some of the early steps taken to tackle this limitation. These 
have also been employed in more genuine non-separable transforms, as seen later, typically relaxing 
one of the aforementioned properties, such as using infinite-support filters [71] , semi- or biorthogonal 
decompositions [59] or complex filter banks [72] . 
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For instance, instead of a two-band filter bank, M-band wavelets [73] with M > 2 provide 
alternatives where symmetry, orthogonality and realness are compatible with finitely supported 
atoms. In this setting, the approximation and the M-band detail spaces are Vj and (Wp) m ^* M 

related through V}_i = Vj @Q^Z\ W™ for a resolution level j. This versatile design provides filters 
that suffer less aliasing artifacts with increased regularity. Their finer subband decomposition is 
also beneficial for detecting orientations in a more subtle fashion than with the {V, H, D} quadrants 
obtained with standard wavelets (Sec. 2.3.2). Yet, more general M-adic MRAs are possible, for 
instance with a rational M = p/q, M > 1 [7H [TSJ [761 E3 EH]- Note that for specific purposes such 
as compression, M-band filter banks with M = 2 J , J e N may be treated like a J-level dyadic tree 
and combined in a hierarchical transform \79\ I80j. Satisfying the MRA axioms is not necessary in 
practice in order to yield high performance results. This is suggested by recent image and video 
coders focusing on "simpler" transforms, closer to ancient Walsh-Hadamard transforms than to 
more involved wavelets [81]. 

Alternatively, the 1-D decomposition on rows and columns of images may be performed in a 
more anisotropic manner, as in [821 [83] . An additional relaxation comes from lifting the critically 
sampled scheme, yielding oversampled, translation-invariant (see Sec. 4.2.3) multiscale wavelets, 
wavelet /cosine packets or frames [671 EH ESI ESI E3 EE] . Multidimensional oversampled filter banks 
in n— D with limited redundancy may be designed as well [891 122 ED 192} [93]. 



3.1.2 Pyramid-related wavelets 

Notably influenced by [9H[95], Unser and Van de Ville propose a slightly redundant transform 
[96] based on a pyramid-like wavelet analysis. This decomposition constitutes a wavelet frame 
with mild redundancy, which is nevertheless not steerable. Subsequently, the same authors propose 
a steerable analysis [97J based on polyharmonic £>-splines [98] and the Maar-like [99] wavelet 
pyramid. Such multiresolution analysis can easily be implemented via filter banks as detailed in 
[97] and the total redundancy of this decomposition is 8/3 (a redundancy of 4/3 is introduced by 
the pyramid structure and the complex nature of the coefficients increases the redundancy by a 
factor of 2). A similar approach based on Riesz-Laplace wavelets is proposed in [100]. The latter 
constructions are related to Hilbert and Riesz transforms. 



3.1.3 Complexifying Discrete Wavelets with Hilbert and Riesz 

Different kinds of complexiflcation are indeed a possible option in order to tackle the problem 
of poor directionality with classical wavelet transforms. The common basic idea leans toward 
analytic wavelets and their combination to improve the 2-D directionality. Behind a generic notion 
of complex wavelets reside different approaches detailed hereafter, which require the definition of 
some basic tools. 

We first introduce the Hilbert transform, termed "complex signal" in [101] and exhaustively 
mapped in [102J. While the 1-D Hilbert transform is unambiguously defined, there exists multidi- 
mensional extensions, often obtained by tensor products, thus leading to approximations. In order 
to increase the directionality property, other multidimensional constructions (discussed in [103] ) 
have also been proposed. 

• The 1-D Hilbert transform % of a signal / is easily expressed in the Fourier domain as 

W}( W ) = -isign( W )/( W ). (6) 
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• The 1-D fractional Hilbert transform %q of / is similarly defined in [104] by 

U e {f}{u) = exp(iyr6>sign(cj))/(cj). 



(7) 



• The 2-D directional Hilbert transform %q of / is one of the 2-D extensions defined in [104] as 
KeifKuii^) = -isign (cos(0)cji + sm(0)u 2 )f(wi,V2)> (8) 

See also [T05] . 

The Hilbert transform was already associated with wavelets for transient detection by Abry 
et al [106] . Others early connections between wavelets and the Hilbert transform are drawn in 
\W7\ \W8\ 1109] . At the end of the 1990's, Kingsbury proposed the dual-tree transform based on 
even and odd filters |110t 1111] . An alternative construction is given by Selesnick [112J. It amounts 
to performing two discrete classical wavelet transforms in parallel, the wavelets generated by the 
trees forming Hilbert pairs. An atom of the corresponding basis (here the diagonal wavelet) and its 
corresponding frequency plane tiling are depicted in Fig. [9| The corresponding dual-tree of wavelet 
coefficients is represented in Fig. |8j which clearly shows the separation of oriented structures with 
different orientations. The resulting oriented wavelet dictionary has a small redundancy and is also 
computationally efficient. The corresponding wavelet is approximately shift invariant, see [113] 
for more details. It is extended to the M-band setting by Chaux et al [ 1T3] and to wavelet 
packets in \115\ I116J . In Fig. [lOj one subband of the wavelet transform (red square in Fig. [7]), 
two subbands (primal+dual) of the dyadic dual-tree transform (red squares in Fig. [8]), as well as 
the corresponding eight subbands (4 primal+4 dual) of the 4-band dual-tree wavelet decomposition 
are depicted. In Fig. |10(d)| the fine oriented textures from the left side of the image are (slightly) 
better separated in some non- horizontal subbands. The wavelet /frequency tiling corresponding to 
the 4-band dual-tree wavelet decomposition are depicted in Fig. [TT] The main advantage of this 
decomposition is that it achieves a directional image analysis with a small redundancy of a factor 
2 (4 for the complex transform). 

Gopinath [117^ 1118] has designed phaselets which is an extension of the dyadic dual-tree wavelet 
transform [110} 1119] . They aim at improving translation invariance with a given redundancy, and 
are built by carefully observing the effects of shifts in a discrete wavelet transform. 2-D phaselets 
are easily obtained by tensor products. 

More recently, the shiftability of the dual-tree transform has been studied by Chaudhury et al 
[104] by introducing the fractional Hilbert transform Q. A 2-D extension has been proposed in 
[120] and the construction of Hilbert transform pairs of wavelet bases can be found in [121] . Note 
that previous works dealing with multidimensional extensions have been first reported for instance 
in [122] and then in \123\ 1124] using the notion of hypercomplex wavelets. 

Numerous extension to multidimensional signals have been proposed, see for instance [12511126] . 
They, for instance, use the Riesz transform 7£, which is defined in the frequency domain as follows: 

£{/} = (£i{/},...,7M/}). (9) 

where 

Vn € {1, • • • , N}, K n {f}(u>) = -iipV/M. (10) 

\\LJ || 

Other recent extensions of multidimensional oriented wavelets are based on the notion of monogenic 
signal/ wavelet [127| [128, 129]. We finally mention that other methods have been developed in order 
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to achieve directional analytic wavelets such as softy space projections [130^ 11314 1132} 1133} 1134] or 
the Daubechies complex wavelets [135^1136^137) . Complex wavelets have also been shown to provide 
robust image similarity measures [138} I139| . 




Figure 8: Dyadic dual-tree wavelet decomposition of the Haar-Riesz Memorial plaque. 
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(a) (b) 

Figure 9: The dyadic dual-tree wavelet, (a) Example of atom (diagonal wavelet), (b) Associated 
frequency partitioning. 



3.2 Non-Separable Directionality 

3.2.1 Non-separable Decomposition Schemes 



In contrast to the separable constructions detailed in Sec. |3.1.1| where n-D representations are 
composed of 1-D transforms applied separately along each dimension (sometimes recombined, as in 
the dual-tree wavelet case or in |140| ). non-separable constructions are directly performed in n-D. 
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(a) 



(c) 




(b) 



(d) 



Figure 10: The original image (a) and the horizontal subband(s) at first resolution level for (b) 
Dyadic wavelet transform, (c) Dyadic dual-tree transform (primal+dual) and (d) M-band dual-tree 
wavelet decomposition (primal+dual) of the Haar-Riesz Memorial plaque. 



(a) (b) 
Figure 11: The M-band dual-tree wavelet, (a) Example of atom, (b) Frequency partinioning. 
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Since the literature on this topic is large, this section is focussed on a limited number of references 
dealing with directional multiscale decompositions. 

These works are often related to non-diagonal subsampling operators, non-rectangular lattices 
(e.g., quincunx grids or integer lattices) [141|. I142j . or non-separable n-D windows [143] 1144] . Com- 
plementary standard references can be found in [701 p. 558 sq.] or [145^ I146[ I147j . Some of these 



constructions are defined using the lifting scheme, see Sec. |4.3| and |5.3| for more details. While di- 
rectional filter banks do not provide a multiscale representation in general, 2-band |148| 11491 1150j or 
even M-band non-redundant directional discrete wavelets [15 lj have been proposed. Non-separable 
schemes are used for instance as building blocks for multiscale geometric decompositions such as: 

• directional filter banks in [152J, and their combination with a Laplacian pyramid in contourlets 
[I531II54] or surfacelets [T55] . 

• (pseudo-) polar fast Fourier transform (FFT) [156] in first generation curvelets described in 



Sec. 3.3.2[ or the loglets in [157] that exhibit a polar separability. 



In order to overcome the limited efficiency of the standard 2-D separable DWT for representing 



non-horizontally or vertically directed edges (see Sec. |2.3.2[ ), several authors have adapted 1-D 
concepts for local edge representation. Reissell [158J develops, for instance, a pseudo-coiflet scheme 
that addresses numerically efficient interpolation for a parametric wavelet representation of curves. 
Moreover, for digital images it would be beneficial to follow contours on more appropriate discrete 
paths (see [159] for an early application) such as discrete lines [16CH 1161] 1162] . While discrete lines 
are adapted to digital ridgelets in [163J, Velisavljevic et al. propose multidirectional anisotropic 
directionlets [164] , based on skewed lattices, with directional vanishing moments along direction 
with rational slopes, still relying on a simple separable implementation. This approach is refined in 
[165] by taking lifting steps of 1-D wavelets along an explicit orientation map defined on a quincunx 
multiresolution sampling grid, and in [166] with a more efficient representation for sharp features. A 
combination of 2-D filter banks and 1-D directional filter bank is devised in [167} 1168] . Similar ideas 
have been recently applied to edge detection in [169] . In [170] . non-adaptive directional wavelet 
frames are constructed with Haar wavelets and a finite collection of "shear" matrices. Krommweh 
also proposes tetrolets, an adaptive variation (akin to digital wedgelets) of Haar-like wavelets on 
compact tetrominoes (geometric shapes composed of four squares, connected orthogonally, see 
[171] ). These last constructions may further sparkle the growing interest of the association of 
multiscale analysis and discrete geometry [172J. 

3.2.2 Steerable Filters 

Steerable filters [173^ 1174^ 1175] were developed in order to achieve more precise feature detec- 
tors adapted to image edge junctions (often termed "X", "T" and "L" junctions). Their construc- 
tion allows one to compute multiscale derivatives at any orientation (steerability) from a linear 
combination of a small number of fixed filters. In [174] . the construction starts from a bidimen- 
sional Gaussian G(x) = exp(— ^||cc|| 2 ) for x = (xi,X2) with associated base (differential) filters 
= ^-G{x) and G^(x) = £-G{x). 

From the properties of the directional derivative, filters "steered" at angle 9 E [0, 2tt) are then 
built from 

G e {x) = cos(e)g°(x) + s\n(6)G*l\x). (11) 
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where cos(#) and sin(0) may be interpreted as interpolators. Since the convolution is linear, the 
resulting steered decomposition arises from a combination of images that underwent Q° or Q 77 / 2 
filters. A larger class of asymmetric oriented filters is proposed in |176| . Their angular parts are 
derived from even and odd functions: 



N 



N 



Mip E [0, 27r), h e ((p) = w n cos(nip) and h (ip) = w n sin(ny), 



(12) 



n=l 



71=1 



which form Hilbert transform pairs (see Sec. 3.1.3), unlike the resulting spatial filters. An angle 9 
rotation is obtained through: 

h e (<p - 9) = k e (6) T f(<p) and h (<p - 9) = k (9) T f(<p), (13) 

where k e (0) and k {9) are interpolating vectors and f((f) is a weighted Fourier vector, namely: 

fc c (0)=[ cos#,sin#, cos(2(9),sin(2#), ••• , cos(7V#), sin(iV#)] T , 
k o {0) = [-sin#,cos#,-sin(2(9),cos(2fl), ••• , - sin(7V#), cos(N9)] T , 
f(ip) = [wi cos ip, w\ sine/?, W2 cos(2(/?), W2 sin(2y?), • • • , wn cos(7V(/?), wjy sm(N(f)] T 

If we set 9 — 9 n — 2tvti/N for 1 ^ n ^ N, filters h e {- — 9) and h Q {- — 9) may be rewritten as a 
linear combination of h e (- — 9 n ) and h Q {- — n ), 1 ^ n ^ N. An example of decomposition with 



four orientations and two scales is represented in Fig. 12, with corresponding projection atoms in 



Fig. 13 Steerable filters may be combined with discrete wavelets to improve their radial properties 

nzaing. 






Figure 12: Steerable pyramid decomposition of the Haar-Riesz Memorial plaque, over two scales, 
with four orientations. 



3.2.3 Directional Wavelets and Frames 



In Sec. 2.2, the two-dimensional Continuous Wavelet Transform (2-D CWT) was defined as a 
straightforward extension of the 1-D CWT using isotropic wavelets. It is however possible to make 
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Figure 13: Example of steerable pyramid atoms with four orientations. 



use of more complicated group actions to drive the CWT parameterization in the plane, such as 
rotations or the similitude group SIM(2), see |147j . 

Consequently, given a mother function E L 2 (M 2 ) that is well localized and oriented, we write 

*l>{b,afi)( X ) = \^-a R e l ( X - h ))i 

where Rq stands for the 2x2 rotation matrix. For a function / E L 2 (R 2 ), the 2-D CWT (non- 
isotropic) is thus 

W f (b,a,0) = (ip {b , a ,e)J)- 

If the wavelet is admissible, i.e., if = (27r) 2 J r2 |^(a;)| 2 /||u;|| 2 d 2 u? < oc, then, the CWT may be 
inverted through 



/(*) 



-4> 



roc r 

7*/ 

Jo Jo 



2tt 

d9 



d 2 b W f (b,a,0) 1> {b ,a,g)(x) 



the equality being valid almost everywhere on Mr. 

The selectivity power of the wavelet, that is, its ability to distinguish two close orientations in 
an image, may be measured in the Fourier domain. Typically, a good directional wavelet is thus 
a function whose Fourier transform is essentially or exactly contained in a cone with apex on the 
origin: the narrower the cone, the more selective the wavelet transform using that wavelet [147J. 

Practically, it is not satisfactory to manipulate a continuum of wavelets parameterized by con- 
tinuous parameters. The question is therefore to know if it is possible to decompose and reconstruct 
an image from a discretized set of parameters, i.e., on the family Q — {V>(&,a,0) b <EV,a E A, E 0} 



with V C M , A C and 9 C [0, 2tt) all discrete (countable) sets. As explained in Sec. 
question amounts to ask when Q is a frame of L 2 (R 2 ). 

Such frames have been built for the Morlet (or Gabor) wavelet |179t 1180 : 



1.3.2 



this 



where coq E 



\\x\\ 2 /2(Jq \<jJqx 



-o;o|| 2 /2 



m 



Fig. 14 



R 2 defines the cone axis and ctq > is related to the cone aperture, as represented 
Notice that approximate quadrature filters exist to accelerate the computation of 
the wavelet coefficients [181J. The Conic (or Cauchy) wavelet, whose spectral support is exactly 
contained into a cone, can also be used in order to define a frame [105J. 

Finally, a multiresolution structure can also be put on the angular dependency of the conic 
wavelets in the frequency domain to define multis elective wavelets |182| . This generates a redun- 
dant basis that may represent jointly a large spectrum of features ranging from highly directional 
ones {e.g., edges) to isotropic elements {e.g., spots, corners) and including intermediate directional 
structures such as textures. 
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(a) 



(b) 



Figure 14: The Morlet Wavelet, (a) Spatial representation (real part), (b) Fourier representation. 
Supporting cone and frequency axes are drawn for illustration. 

3.3 Directionality in Anisotropic Scaling 
3.3.1 Ridgelets 

Ridgelets \183\ 1184] and wavelet X-ray transforms [185] appear as a combination of a 1-D 
wavelet transform and the Radon transform [186] . They are designed for efficient representation of 
discontinuities over straight lines. A bivariate ridgelet transform is constant along parameterized 
lines x\ cos(#) + X2 sin(0) = b and defined for a > 0, b E K and 9 E [0, 27r), by 



= (xi,x 2 ) E M 2 , ^(b,a,0)( x ) = a 1 ^ 2 '0((^icos(6 ) ) + x 2 sin(60 - b)/a) 
Ridgelet coefficients for the image / are given by 

n f (b,a,0) = J \p M) (x)f(x)d 2 x 

= JlR f (0,t) a- 1/2 t/j((t-b)/a) di, 

where !EH/(#,£) represents the Radon transform of / defined by: 

9\f(6,t) = J J f(xi,X2)S(xicos(0)+X2sm(0) — t) dxidx 2 , 



(14) 



(15) 



(16) 



with 5 denoting the Dirac distribution. The ridgelet transform may be interpreted as a 1-D wavelet 
transform of Radon slices where the angle 9 is constant and t varies. Several implementations 
and variations exist in order to overcome the issues raised by the Radon transform discretization, 
such as the finite ridgelet transform |187| . the approximate digital ridgelet transform [188] or the 
discrete analytical ridgelet transform [189] . Their multiscale implementation [16] is the basis for the 
first generation curvelets described in Sec. 3.3.2 A ridgelet decomposition 7 [190] of the Haar-Riesz 



Memorial plaque is given in Fig. 15, with a typical atom along with a synthetic description of its 
implementation in Fig. [TBI 



BeamLab toolbox: http://www-stat.stanford.edu/~beamlab/ 
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Figure 15: Ridgelet decomposition (square root scale) of the Haar-Riesz Memorial plaque. 



3.3.2 Curvelets 

The curvelet representation, introduced by Candes and Donoho [TH I191j , improves the approx- 
imation of cartoon images with C 2 edges with respect to wavelets. We review here the second 
generation of curvelets, as introduced in [191 j . 

Continuous Curvelet Transform A curvelet atom, with scale s, orientation 9 E [0, 7r), position 
y E [0, l] 2 is defined as 

tl>8,yA*) = MR* 1 (*-v)) (17) 

where i/j s (x) w s -3 / 4 ^(5 _1 / 2 xi, s~ 1 X2) is approximately a parabolic stretch of a curvelet function 
i/j with vanishing moments in the vertical direction. At scale s, a curvelet atom is thus a needle 
oriented in the direction 9 whose envelope is a specified ridge of effective length s 1 / 2 and width s, 
and which displays an oscillatory behavior transverse to the ridge. 

A curvelet atom thus benefits from a parabolic scaling property width = length 2 that is a major 
departure from oriented wavelets. Fig. [17] presents an example of a curvelet atom, together with 
its Fourier transform, for the second generation of curvelets. The resulting curvelet Fourier tiling 
resembles that of the Cortex transform [192J. 

The continuous curvelet transform computes the set of inner products (^ s ,y,o(')i f) f° r an possi- 
ble (s, 2/, 9). A careful design of i(j s (191 j enables a conservation of energy and a simple reconstruction 
formula. The decay of the curvelet transform as s decreases allows one to detect the position and 
orientation of contours [ 193j . 

Curvelet Frame The continuous curvelet representation is sampled in order to obtain a curvelet 
frame B — {^ m } m , |191| . see also |194| for the description of a complex curvelet tight frame. 
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(a) 




Radon Transform 



Ridgelet 



Frequency 



(b) 

Figure 16: The Ridgelet transform, (a) Example of atoms, (b) Synthetic implementation descrip- 
tion. 



A curvelet atom, with scale 2 J , orientation 0£ E [0,7r), position x n E [0, l] 2 is defined from the 



continuous atom (17) 



where the sampling locations are 

^=^7r2^"/2J-l e [ 0j7r ) and x n = R 0e (2 j / 2 n u 2% 2 ) E [0,1] 2 . 

The curvelet parameters are sampled using an increasing number of orientations at finer scales. 
This sampling is the key ingredient to ensure the tight frame property |191| , which provides a 
simple reconstruction formula. 

A fast discrete curvelet transform computes the set of inner products {(^ m , f)}m m 0(N \og(N)) 
operations for an image with N pixels, see [195J. The coronae and rotations of the continuous set- 
tings are replaced by their discrete Cartesian counterparts, i.e. concentric squares and shears. 
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Figure 17: Left: Example of a curvelet ^(cc,«s,y, 
is a wedge. 



). Right: the frequency support of ^(cj, s, j/, 0) 



Figure [18] shows an example of curvelets decomposition . 

Candes and Donoho prove [196J that the curvelet non-linear approximation /m = St(/ 5 £>), 
where Ht is defined in Q, ensures an approximation error decay ||/ — /m|| 2 = 0(M~ 2 log 3 (M)) 
for a C 2 regular image outside C 2 regular edge curves. This is a significant improvement over the 
0(M~ 1 ) error decay of a wavelet approximation described in Sec. 



2.3.2 



and is achieved with a 

fast 0(N\og(N)) algorithm for discrete images. This asymptotic error decay is optimal (up to 
logarithmic factor) for the class of images that are C 2 regular outside C 2 regular edge curves, see 
|196j . Monogenic curvelets are proposed in |197j to obtain additional advantages over monogenic 
wavelets, described in Section [3.1.3 



Shearlet atoms [198 ; 199J are built similarly to curvelets, but they replace, in their continuous 
formulation, rotation and anisotropic stretch with anisotropic shears. The discrete shearlet trans- 
form |200t 1201] is thus implemented similarly to the discret curvelet transform [195] using discrete 
shears 9 . It provides the same approximation properties as curvelets, albeit with a different direc- 
tional sensitivity (e.g., the number of orientations doubles at each scale). Recently a type-I ripplet 
transform [202J has been proposed as an extension to curvelets with alternative scaling laws. 



3.3.3 Contourlets 



Contourlets |153j are sometimes considered a low-redundancy discrete approximation of curvelets. 
Actually, they are designed in the spatial domain (instead of the frequency plane), aiming at a close- 
to-critical directional representation. Their construction is based on a Laplacian Pyramid [13] (see 
Fig. [5]). The low-pass part of the pyramid is further decomposed with a biorthogoal 9/7 DWT. 
Each difference image obtained from the pyramid is subject to directional filter bank (see Sec. |3.2[ ) 
(initially from [141] . [203J proposes a simpler implementation based only on a quincunx structure). 
A contourlet decomposition is illustrated 10 in Fig. 19, The resulting frequency plane tiling is rep- 
resented in Fig. |20(c) The contourlet inherits its redundancy of 4/3 from the pyramidal scheme. 
Its approximation rate is similar to that of curvelets (Sec. 3.3.2). At one end of the redundancy 



8 The Curvelab toolbox has been u sed, see |http : // www . curvelet . org/[ 

9 An implementation is available at http : / / www . shearlab . org 

10 The contourlet toolbox has been used, see Ihttp : //www . if pTillinois . edu/~minhdo/sof tware/ 
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Figure 18: Curvelet decomposition of the Haar-Riesz Memorial plaque. The layout of the coeffi- 
cients follows the frequency localization of curvelet atoms. 
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spectrum, |204j proposes a critically sampled version. At the other end, the constraints thus laid on 
the basis functions (Figs. 20(a)p0(b) ) are relaxed by the design of a more redundant [ 154| version, 
based on non-subsampled (Sec. 3.1.1 )pyramid and directional filters. 
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Figure 19: Contourlet decomposition of the Haar-Riesz Memorial plaque. 




3.3.4 Frames for Oscillating Textures. 

While curvelets, contourlet s and shearlets are optimized for the processing of edges, they are 
not tailored for the processing of oscillating textures, because of their poor frequency localization. 
Generic oscillating patterns can be captured using a local Fourier analysis on a regular segmentation 
of the image in squares. This corresponds to an expansion in a Gabor frame, see for instance [33] . 
The spatial segmentation can be optimized using a decomposition in a best cosine packet dictionary 
as described in Section 14.21 



Wavelet packets, detailed in Section 4.2, have been used to process and compress oscillating 



textures such as fingerprints. Brushlets [205J, introduced by Meyer and Coifman, improve the 
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frequency localization of wavelet packets. 

Wave atoms [206J better capture geometric textures using an anisotropic scaling 11 . The wave- 
length of wave-atom oscillations is proportional to the square of their diameter. This scaling allows 
a thresholding in a wave atom frame to optimally approximate textures obtained by a smooth 
warping of a sinusoidal profile, see [206J. 



4 Redundancy and Adaptivity 

Highly redundant representations allow us to improve the representation of complicated images 
with edges and textures. However, as described hereafter, computing efficient image representations 
in such dictionaries sometimes requires approximations. 



4.1 Pursuits in Redundant Dictionaries 

An approximation /m of an image / with M atoms from a highly redundant dictionary B — 
{^mj ' 1 < j < P} is written 

f M = = ^2 a j with ||a|| = # { j '• clj ^ 0} < M. 

3 

Computing the M-sparse coefficients a that produce the smallest error ||/ — /m\\ in a generic 
dictionary is NP-hard [207] . Furthermore, the M-terms approximation Jm — Ht(/^B) computed 
by thresholding ^ might be quite far from the best M-terms approximation. One thus has to use 
approximate schemes in order to compute an efficient approximation in a reasonable time. 



4.1.1 Matching Pursuits 

Matching pursuit [208] computes /m from /m-i by choosing the atom ^ m that minimizes the 
correlation K^m, / — /m-i)|- Orthogonal matching pursuit [33 } 1209] further reduces the approxi- 
mation error by projecting / on the M chosen atoms to compute /m- 

Under restrictive conditions on the dictionary £>, these greedy algorithms compute an approx- 
imation /m that is close to the best M-term approximation, see for instance [2101 211J. These 
conditions typically require the correlation |(^ m , Vw)| to be small for m ^ m', which is not 
applicable to highly redundant dictionaries typically used in image processing. 



4.1.2 Basis Pursuit 



A sparse approximation is obtained by convexifying the fi N pseudo norm, and solving the 
following basis pursuit denoising convex problem [212] 



f M = \Jj a = ^aj^rnj where a e argmin ±\\f - ^ a^^. || 2 + ^\\a\\\, 



(18) 



where ^ > is adapted so that ||a||o = M. This problem (18) is minimized, for instance, using 
iterative thresholding methods \213\ 1214] . Algorithmic solutions to its generalized form as sums of 



L See 



http : //www . waveatom . org 
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convex functions (a common formulation to many data processing problems) may be solved with 
great flexibility in the framework of proximity operators [215J. 

Similarly to matching pursuit algorithms, this i]^ approximation can be shown to be close to 
the best M-term approximation if the atoms of B are not too correlated, see for instance |216}I211| . 

4.1.3 Pursuits in Parametric Dictionaries 

Parametric dictionaries are obtained from basic operations (like rotation, translation, dilation, 
shearing, modulation, etc.) applied to a continuous mother function. Even if such dictionaries 
also define redundant bases similar to those introduced earlier, they deserve a separate description 
since their parametric nature provides them with some particular properties. They are generally 
created to provide a very rich and dense family of functions built from the geometrical features of 
the analyzed image. They have applications in image and video coding [217J, multi-modal signal 
analysis (e.g., video plus audio) |218j . and also for signal decomposition on non-Euclidean spaces 

mm. 

Formally, given a set of S transformations . for 1 ^ i ^ S parameterized by rrii E C , 
the parametric dictionary is related to a certain discretization of A d C A = Ai x • • • x A#, i.e., 

B = {Vm(tf) = [T^ ■ ■■T^]{x) € L 2 (M 2 ) : m = (mi, • • • , m s ) € A d }. 



The directional wavelets described in Sec. |3.2.3 and the subsequent frames built from them 



are actually an example of parametric dictionaries with the translations Z^Z^ = T£ t£ , the 
rotation T^ 3 = Rq and the dilation ZJ^ 4 = D a operations. For these wavelets, the decompo- 
sition/reconstruction methods are relatively easy to formulate, due to the continuous inversion 
formula or using the frame condition. 

However, checking the frame condition may sometimes become tedious. In addition, more 
transformations of the mother function may be added in order to enlarge the family of functions, 
further worsening the frame bounds. 



Best of A 



Best of A C 




Figure 21: Explanation of the optimization in A starting from a point in A d . 



Fortunately, as described in Sec. |4.1| it is still possible to find good description of images in very 
general family of functions. Most of the time, since the Parametric Dictionaries are much larger 
than other dictionaries of controlled redundancy, the (Orthogonal) Matching Pursuit decomposition 



(Sec. 4.1.1) is used to find a sparse representation of signals. 
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(a) 



(b) 



w P L 



(c) 



Figure 22: (a) Original image, (b) Reconstruction with 300 atoms for a rich parametric dictionary 
containing 5x5 anisotropic scales, 8 orientations, and N translations. PSNR : 26.63 dB (CT: 4634s). 
(c) Optimized Reconstruction at 300 atoms starting from a dictionary with only 3x3 scales, 4 
directions and N translations, PSNR: 26.68 dB (CT: 949s). 



Interestingly, thanks to the parametric nature of £>, the dictionary discretization can be refined 
during the Matching Pursuit iterations. Indeed, since B is the discretization of the continuous 
manifold M = {ipm : m E A} C L 2 (M 2 ) generated by all the transformations of ^, at each 
iteration of MP in the decomposition of a signal / E L 2 (R 2 ) the refinement is performed as follows. 
As illustrated on Fig. [21} given the best atom ^ m found in £>, a gradient ascent respecting the 
(Riemannian) geometry of A4 is run on A to maximize the correlation S(m f ) = |(Vw 3 RJ)\ between 
the current MP residual RJ = f — f n at step n and the atom ^W- A new parameter m* is then 
used instead of m in the signal representation and the next iteration is realized on the residual 

R*}-(^rn*,R n f )^rn* [220j . 

Fig. 22 presents the result of such an improvement for two different decompositions of the 
Barbara image (with N = 128 2 pixels) with similar qualities (expressed using the Peak Signal-to- 
Noise Ratio - PSNR). The first one (Fig. |22(bJ| ) is obtained by a rich parametric dictionary defined 
by anisotropic dilations, rotations, and translations of a 2-D second order directional derivative of 
a Gaussian. The second decomposition uses a poorer dictionary with the same parameterization 
and mother function but with a manifold optimization on the atom parameters. The interest of 
the latter method is to provide a similar quality for a smaller Computational Time (CT). 

4.1.4 Processing with Highly Redundant Dictionaries 

Compression with Sparse Expansions Dictionaries with oriented atoms have proven to be 
successful for improving the JPEG 2000 compression standard at low bit rates |221( 1222] . The 
approximation of the image is computed using the matching pursuit algorithm. Matching pursuit 



in Gabor dictionaries, i.e., dictionaries made of Gabor wavelets (Sec. 3.2.3), have been used for 
coding the motion residual in video compression schemes [223J. 



Inverse Problem Regularization Data acquisition devices usually only acquire S noisy low 
resolution measurements y = <&/o + w E of a high resolution image fo E M, N of N ^> S pixels. 
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(a) 



(b) 



(c) 



Figure 23: Example of deconvolution using t x N regularization in a frame of translation invariant 
wavelets, (a) Original /q. (b) Observation y = $/o + (c) Deconvolution /. 



The linear operator $ models the acquisition and might include some blurring and sub-sampling 
of the high resolution data. 

Recovering a good approximation / E M. N of f$ from these measurements y corresponds to 
solving a difficult ill-posed inverse problem, that requires the use of efficient priors to model the 
regularity of the image. Early priors include the Sobolev prior that enforces smoothness of the 
image, and the non-linear total variation |224j that can produce sharper edges. 

More recently, l x N sparse priors in redundant dictionaries B have been proved to be efficient in 
order to solve several ill-posed problems, see for instance [33] and references therein. In this setting, 
one computes the coefficients a of / = \I/a = Y2j a j 1 Pmj hi a frame B of P atoms by solving a l x N 
augmented Lagrangian form 



a E argmin ^\\y 



$*a|| 2 + n\\a\\ 



where \I/a = ajt(j r 



(19) 



where /i should be adapted to the noise level ||n;|| that is supposed to be known. This minimization 
problem corresponds to computing the basis pursuit approximation (18) of the measurements y in 
the highly redundant dictionary {<3>^ m . : 1 ^ j ^ P} of W s . It can thus be solved using the same 
algorithms. 



Figure 23 shows the use of this sparse regularization method when solving a deconvolution 
problem. In this application, the operator is a convolution <&/ = f*G a with a Gaussian kernel G a 
as defined in Sec. |1.3.1| The redundant dictionary B is a translation invariant wavelet frame. 



4.1.5 Source Separation 

Sparse representations can be used to separate sources that are known to be sparse in different 
dictionaries. This corresponds to the morphological component analysis (MCA) of Starck et al 
[225] . In its simplest setting, it can be used to separate a single noisy image y into a sum y = 
Jg + It + w of a cartoon-like component fc (or geometric component), a texture component 
and residual noise w. One can use a dictionary B = Bq U Bt union of wavelets (Bq) and local 
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(a) (b) (c) 



Figure 24: Example of cartoon+texture decomposition using the MCA algorithm, (a) Original y. 
(b) Geometry layer /q. (c) Texture layer fa. 



cosine (£>t) 5 an d compute a sparse approximation / of y 

f = *a = ^c a G + ^tclt 



(20) 



where a = [a<^; aj-] is the solution of the fl N basis pursuit ( fl8| ) applied to y. The separation, obtained 
using fc = ^ G a G and /j- = ^t^Ti is illustrated in Fig. [24j 

The modeling of natural images as a sum of a cartoon layer and an oscillating texture layer has 



been initiated by Y. Meyer in his book [12] • Beside sparsity-based approaches such as (20), other 
variational methods have been proposed, see for instance the work of J.-F. Aujol et al [T3] . 



4.2 Tree-structured Best Basis Representations 

Pursuit algorithms are quite slow and face difficulties in order to compute provably efficient 
approximations when the dictionary is too redundant. In order to avoid these bottlenecks, one 
needs to consider more structured representations, that allow one to use fast and provably efficient 
approximation strategies. The structuring of the representation can be implemented by computing 
an adapted basis B x parameterized by a geometric parameter A that captures the local direction 
of edges or textures. This section details best basis schemes: they introduce the desired adaptivity 
together with fast algorithms employing the hierarchical structure of parameters A. 



4.2.1 Quadtree-based Dictionaries 

A dictionary of orthonormal bases is a set T)\ — agA of orthonormal bases B x — {^Prn}m 
of M N i where N is the number of pixels in the image. Instead of using an a priori fixed basis such 
as the wavelet or Fourier basis, one chooses a parameter A* £ A adapted to the structure of the 
image to process and then uses the optimized basis B x . 

In order to enable the fast optimization of a parameter A* adapted to a given signal or image 
/ to process, each A £ A is constrained to be a quadtree. The quadtree A that parameterizes a 
basis B x defines a dyadic segmentation of the square [0, l] 2 = \J(ji)eL(\) where L(A) are the 



leaves of the trees, as shown on Fig. 25, Each square Sj^ is recursively split into four sub-squares 
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Sj+i^i+k for fc = 0, • • • , 3. In order to enrich the representation parameterized by a quadtree, we 
attach to each leave of the tree a geometric token, and denote as r the number of tokens. A token 
indicates the direction of the image geometry in a square of the segmentation. 

4.2.2 Best Basis Selection 

Given a number M of coefficients, the best basis £> A * E V\ adapted to / E M. N minimizes the 
best M-terms approximation error. This can be equivalently obtained by minimizing a penalized 
Lagrangian that weights the approximation error with the number of coefficients 

A* E argmin £(/, B\ T) = \\f - f^f + M A T 2 , (21) 

where fj^j is the best M A -term approximation in B x computed by thresholding at T > 

f&=H T (f,B\T)= J2 and M A = # {m : |(^, /)| > t] , (22) 

K^,/)|>r 

since A is orthonormal. This Lagrangian can be re-written as a sum over each coefficient in the 
basis 

C(f, B X ,T) = Y, max(| (^i, /) | 2 , T 2 ). (23) 

m 

This kind of Lagrangian can be efficiently optimized using a dynamic search algorithm, originally 
presented by Coifman et al [226J, which is a particular instance of the Classification and Regression 
Tree (CART) algorithm of Breiman et al |227j as explained by Donoho [228] . It is possible to 
consider other criteria for best basis selection, such as for instance the entropy of the coefficients. 
This leads different Lagrangians that can be minimized with the same method [226J. 

The complexity of the algorithm is proportional to the complexity of computing the whole 
set of inner products {(^m? /) : A E A} in the dictionary. For several dictionaries, such as those 
considered in this section, a fast algorithm performs this computation in 0(P) operations, where 
P is the total number of atoms in V\. For tree-structured dictionaries, this complexity is thus 
0(rAnog 2 (AQ), where r is the number of tokens associated to each leaf of the tree. This is much 
smaller than the total number of basis £> A in that grows exponentially with N. 
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4.2.3 Wavelet and Cosine Packets 



A basis B x with oscillating atoms is denned using a separable cosine basis over each square of 
the dyadic segmentation. In this case no geometry is used, the oscillation of the atoms does not 
follow the geometry of the image, and r — 1. An approximation in an adapted cosine basis B x 
allows one to capture the spatial variations of a texture [33] . 

A wavelet packet basis B x defines a dyadic subdivision of the 2-D frequency domain [229] . The 
projection of an image on the atoms of B x is computed through a pyramidal decomposition that 
generalizes the orthogonal wavelet transform, adding flexibility to overcome its dyadic frequency de- 
composition. Uniform dyadic wavelet packet decompositions generate a subset of M-band wavelets 
with equal-span frequency subbands obtained from J decomposition levels, with M — 2 J . In order 
to adapt to the specific frequency content of the image, the resulting tree is parsed through a best 



basis selection procedure |226| . reminiscent of the subdivision in Fig. 25 



This construction is generalized by considering non- stationary wavelet packets [230J , that apply 
different quadrature mirror filters at each scale of the tree. A dynamic programming algorithm 
detailed in [231J computes an adapted non-stationary basis. 



4.2.4 Adaptive Approximation 

Wedgelets A geometric approximation is obtained by considering for each node of the dyadic 
segmentation a collection of r different low-dimensional discontinuous approximation spaces [232J. 
For each node of the quadtree, a token indicates the local direction and position of the edge. The 
low-dimensional approximation spaces are piecewise polynomials over each of the two wedges. 

The wedgelets introduced by Donoho [232J rely on piecewise constant approximation. This 
scheme is efficient when approximating a piecewise constant image / whose edges are C 2 curves. 
For such cartoon images, the approximation error decays like ||/ — /m|| 2 — M -2 , see |232^ [2T]. It is 
also possible to consider approximation spaces with higher-order polynomials in order to capture 
arbitrary cartoon images [233] . see also [234j for a related construction. The computation of the 
low-dimensional projection can be significantly accelerated, see [235j . 

The piecewise constant model for images being relatively simplistic, wedgelets have been up- 
graded to platelets [236J and surflets [237J. They aim at improving the management of smooth 
intensity variations, since they rely on planar or even smoother approximation on dyadic square or 
wedge based grids. 



Bandlets For coding, orthogonal expansions are preferred over low-dimensional approximations 
as considered by wedgelets. Switching to non-linear approximation in bases also better handles 
directional textures that do not correspond to a fixed low-dimensional space parameterized by a 
wedge. 

The bandlet bases dictionary is introduced by Le Pennec and Mallat [238J. Bandlets perform 
an efficient adaptive approximation of images with geometric singularities. An anisotropic basis 
with a preferred orientation is defined over each square of the dyadic segmentation. Fig. [26] (a) 
shows an example of bandlet atom. The orientation is parameterized with the token stored in the 
leaf of tree. Keeping only a few bandlet coefficients and setting the others to zero performs an 
approximation of the original image that follows the local orientation indicated by the token. 



33 



(a) (b) 

Figure 26: Example of a bandlet atom, (a) Atom in the spatial domain, (b) Wavelet-bandelet 
atom. 



Adaptive Approximation over the Wavelet Domain Applying such an adaptive geometric 
approximation directly on the image leads to unpleasant visual artifacts. In order to overcome this 
issue, one applies a tree-structured approximation or a best basis computation on the discrete set 
of wavelet coefficients. The wedgeprint of Wakin et al |239| uses a vector quantization to extend 
the wedgelet scheme to the wavelet domain. The second generation bandlets of Peyre and Mallat 
[240] use an adaptive bandlet basis for each scale of the wavelet transform. All these methods 
benefit from the same approximation error decay as their single scale predecessors, but work better 
in practice. 

Fig. 26 shows how a bandlet atom (a) is mapped to a wavelet-bandlet atom (b). Decomposing 
an image over a bandlet basis composed of atoms of type (b) is equivalent to applying first a wavelet 
transform, and then decomposing the wavelet coefficients over atoms of type (a). 

Another adaptive approximation relying on the processing of the wavelet domain is the easy 
path wavelet transform (EPWT) [241] . It provides a hybrid and adaptive approach exploiting the 
local correlations of images along path vectors through index subsets in the Wavelet domain. 



4.2.5 Adaptive Tree-structured Processing 

For compression and denoising applications, one computes the best basis £> A * adapted to the 
image / to compress or denoise by minimizing the corresponding Lagrangian ( [23] ) . The coefficients 
(^m? /) are then binary coded (for compression) or thresholded (for denoising). The resulting 
improvement of the best basis approximation error over wavelets translates into improvement in 
the rate distortion (for compression) or average risk (for denoising) of the best basis method, see 
for instance [2391 [M)] . 

One can also use best bases to recover an image from noisy low-dimensional measurements 
y = <$>f + w where $ is an ill-conditioned linear mapping. For some problems such as inpainting, 
small missing regions or light blur removal, the best basis A can be estimated directly from the 
observation y. 

An example of inverse problem where sparsity in a best basis significantly improves over sparsity 
in a fixed basis is compressed sensing. Compressed sensing is a new data sampling strategy, where 
the measurement operator <3> of size P x N is generally the realization of some random matrix 
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Figure 27: (a,a') original image ; (b) compressed sensing reconstruction using a translation invariant 
wavelet frame (PSNR=37.1dB) ; (c) reconstruction using a best bandlet basis (PSNR=39.3dB). 
(b') wavelet frame, PSNR=22.1dB, (c') bandlet basis, PSNR=23.9dB. 



ensemble. The sampling operations y = <&/ + w E M p allows one to acquire a high resolution signal 
/ E M. N directly in a compressed format of P < N measurements. Compressed sensing theory 
ensures that if the number of measurements P is large enough with respect to the sparsity K of the 
signal / in a basis £>, typically, P — 0{K log N/K) for Gaussian random matrix <3>, one recovers 
a good approximation of the signal using a t x N sparse regularization as in (19). It can be shown 
that the quality of the reconstruction depends both on the sensing noise power \\w\\ and on the 
"compressibility" of /, that is, its deviation from the strictly sparse case. We refer to the review 
paper of Candes |242j and the references therein for more details. Fig. 27 shows a comparison of 
compressed sensing recovery from P — N / 6 measurements using a redundant frame B of translation 
invariant wavelets, and a best bandlet basis. In this last result, it is necessary to use an iterative 
algorithm that progressively improves the quality of the estimated geometry, see [243] . As explained 
in this last reference, the same technique can be used for inpainting large holes in images. 



4.2.6 Adaptive Segmentations and Triangulations 

In order to enhance the quality of the representation, it is possible to consider tree-structured 
segmentations [0, l] 2 = U/3eA P °f the image where the boundaries of the sub-domains E A are not 
restricted to be axis-aligned. The advantage is that such an adaptive segmentation defines regions 
(3 E A with arbitrary complicated boundaries. Unfortunately, the combinatorial explosion of the 
set of all possible A forbids the search for an optimal segmentation with a fast algorithm. One has 
thus to use a greedy scheme that selects at each step a split to reduce the approximation error. 



Recursive Splitting and Approximation Spaces A greedy scheme computes an embedded 
segmentation A = {Aj}j, where A J+ i C Xj is obtained by splitting a region /3 E Xj. The full 
segmentation A can thus be represented and coded using a binary tree. This defines multiresolution 
spaces V\ j+1 C V\ j where V\ j is composed, for instance, of piecewise polynomial functions on each 
region (3 E Xj. 
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It is possible to compute a single-scale orthogonal projection /m = Pv x .(f) °f an i ma g e / on 
a fixed resolution space V\ j in order to perform image approximation or compression. It is also 
possible to define a detail space V\ j+1 = V\ j © W\ j . A wavelet basis B x can be built by considering 
a basis for each W\ r A non-linear thresholding approximation Jm — ^t(/ 5 £> a ,T) provides an 
additional degree of adaptivity and reduces the approximation error ||/ — /m\\- Wavelet bases on 
adaptive segmentations also enable a progressive coding of the coefficients by decaying T, which is 
important for image compression applications. 



Adaptive Segmentation A popular splitting rule is the binary space tiling, that splits a region 
(3 E Xj according to a straight line, see for instance |244j . 

Other popular approaches restrict the regions f3 E Xj to triangles, so that Xj is a triangulation 
of the domain [0, l] 2 . It is possible to refine the triangulation by adding new vertices, or on the 
contrary to remove vertices to go from Aj+i to Xj. These vertex-based schemes do not satisfy 
Aj+i C Aj, so one cannot build a wavelet basis using such triangulations. These vertex refinement 
methods generate a single scale approximation Py x (/) and lead to efficient image coders, see for 
instance |245j . 

To generate embedded approximation spaces Aj+i C Aj, one needs to split the triangles (3 E Xj. 
Regular split of orthogonal triangles leads to isotropic adaptive triangulations [246]. Splitting 
triangles according to a well chosen median leads to anisotropic triangulations that exhibit optimal 
aspect ratio for smooth images, see [247J. More complicated, non-linear coding schemes are possible, 
for instance using normal meshes [248J, that treat an image as an height field. 



4.3 Lifting Representations 

To enhance the wavelet representation, the wavelet filters can be adapted to the image content. 
The lifting scheme, popularized by Sweldens |249| and latent in earlier works [2501 12511 1252] , is an 
unifying framework to design adaptive biorthogonal wavelets, through the use of spatially varying 
local interpolations. While it can typically reduce the computation of the wavelet transform by a 
factor of about two in 1-D, it also guarantees perfect reconstruction for arbitrary filters, and can 
be used (Sec. 5.3) on non-translation invariant grids to build wavelets on surfaces, see Sec. [51 



4.3.1 Lifting Scheme 

At each scale j, the scaling coefficients ctj-i are evenly split into two groups a° and d°-. The 
wavelet coefficients dj and the coarse scale coefficients aj are obtained by applying linear operators 
P^ 3 and Uj 3 parameterized by Aj 

dj = d° 3 - p] 3 a° 3 and aj = a] + \j] 3 dj . (24) 

The resulting lifted wavelet coefficients {dj}j are thresholded or quantized to achieve denoising or 
compression. These two lifting or ladder steps are easily inverted by reverting the order of the 
operations. The predictor P^ 3 interpolates the sub-sampled values a° in order to reduce the am- 
plitude of the wavelet coefficients dj, while the update mapping Uj 3 stabilizes the transform by 
maintaining certain quantities such as the mean of the scaling coefficients. By applying sequen- 
tially several predict and one update operators, one can recover arbitrary biorthogonal wavelets on 
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(a) 



(b) 



Figure 28: (a) Predict and update lifting steps (b) MaxMin lifting of the Haar-Riesz Memorial 
plaque. 

uniform 1-D grid [253], speeding up the wavelet decomposition algorithm by a factor of about two 
in 1-D. The lifting structure in Fig. |28(a)| corresponds to the 5/3 lifted wavelet. Such structures 
may furthermore adapt to non-linear filters and morphological operations [254 ; 255J. An example 12 
of lifting based quincunx scheme example from \256\ I257j is displayed in Fig. |28(b)[ 

4.3.2 Adaptive Predictions 

It is possible to design the set of parameter A = {Xj}j to adapt the transform to the geometry 
of the image. We call Xj an association field, since it typically links a coefficient of a° to a few 
neighboring coefficients in d°. Each association is optimized to reduce, as much as possible, the 
magnitude of wavelet coefficients dj, and should thus follow the geometric structures in the image. 
One can compute these associations to reduce the length of the wavelet filter near the edges, using 
the information from the coarser scale [258] . Locally adaptive schemes have proven efficient in 
stereo and video coding [2591 EH l26Tll262j . 

Such schemes are related to adaptive non-linear subdivision [263]. To further reduce the distor- 
tion of geometric images, the orientations of the association fields {Xj}j can be optimized though 
the scales. Because of the lack of structure of the set of bases £> A , computing the field Xj that 
produces the best non-linear approximation is intractable. These flows are thus usually computed 
using heuristics to detect the local orientation of edges, see for instance \264\ 1265] I165( 1266] . These 
adaptive lifting schemes are extended to perform adaptive video transforms where the lifting steps 
operate in time by following the optical flow Aj, see for instance |267i 268J. 

4.3.3 Grouplets 

A difficulty with lifted transforms is that they do not guarantee the orthogonality of the resulting 
wavelet frame. The stability of the transform thus tends to degrade for complicated association 

12 LISQ toolbox: http: //www.mathworks . com/matlabcentral/f ileexchange/ 13507 
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fields {Aj}j. The grouplet transform, introduced by Mallat [269J, also makes use of association fields, 
but it replaces the lifting computation of wavelet coefficients by an extended Haar transform, where 
coefficients in d°- are processed in sequential order to maintain orthogonality. 

Grouplets defined over each scale of the wavelet transform have been used to perform image 
denoising, super-resolution [269J and inpainting |270j by solving a l x N regularization similar to (19). 



Grouplets can also be used to solve computer graphics problems such as texture synthesis. Classical 
approaches to texture synthesis use statistical models over a fixed representation such as a wavelet 
basis, see for instance |271j I272j . Building similar statistical models over a grouplet basis |270j 
allows one to better synthesize the geometry of some textures, and gives results similar to state 
of the art computer graphics approaches such as texture quilting [273J. Furthermore, the explicit 
parameterization of the geometry though the association fields A allows the user to modify this 
geometry and synthesize dynamic textures. A comparison of these different approaches on one 



texture synthesis example is given in Fig. 29 




(a) 



(b) 



(d) 



Figure 29: Example of texture synthesis by statistical modeling of grouplet coefficients, (a) Exem- 
plar, (b) Wavelet [272]. (c) Quilting [273]. (d) Grouplets [270]. 



5 Transformations on Non-Euclidean Geometries 

In this section we describe how the concepts of frequency, scale and even directionality have 
been extended to the processing of data on non-euclidean geometries like the sphere and other 
manifolds. 

5.1 Data Processing on the Sphere 

The unit sphere S 2 — {x E R 3 : ||x|| = 1} C K 3 is one of the most natural non-Euclidean spaces. 
Very early, possibly due to influences for astronomy and geosciences, many data processing tech- 
niques have been developed for this surface. Many filtering, multiscale, directional and hierarchical 
methods have been designed, either in the spherical frequency domain induced by the spherical 
harmonics basis — often following the spirit of some Euclidean techniques exposed in the previ- 
ous sections — or on the sphere itself thanks to some geometrical tools such as the stereographic 
dilation or the lifting schemes for wavelet analysis. 
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5.1.1 Filtering 

As for the plane, filtering operations may be defined on S 2 . Given the common two-angle 
spherical parameterization £ = (0, (p) E S 2 with the co-latitude 9 E [0, tt] and longitude ip E [0, 27r), 
this operation is realized through spherical convolution evaluated on SO (3) (the group of rotations 
in R 3 ). For a function / E h 2 (S 2 ) = {g : \\g\\% = f s2 \g\ 2 < oc} and a filter h E h 2 (S 2 ), the 
convolution is 

(f*h)( P ) = [ MpO/(0 dMO, 

Js 2 

where p E SO(3) is a rotation (driven by three angles) applied to the point £ E S 2 and d/i(0 = 
sin#d#d(/?. For an axisymmetric filter, i.e., if h(£) = h(0), the convolution reduces to (/ * h)(£ f ) = 
J S 2 h(£' • 0/(0 d/i(0? where £' • £ is the common 3-D scalar product between £' and £ seen as 
unit vectors. 

5.1.2 Fourier Transform 

The Fourier transform of a function / E JL 2 (S 2 ) is defined by 

hm = (Y £m , f)= f y; m (0 /(0 d//(0 , /(I) = E /*» y ^(0 

with respect to orthonormal basis of spherical harmonics y = {Y£ m (^) : £ ^ 0, |m| ^ £}, i.e., the 
eigenvectors of the spherical Laplacian [274] . 

The frequency content of / is thus represented by the value of fi m on the order £ E N, which 
basically counts the number of oscillations on the latitudes, and the moment m E {—£,- — , £} 
counting longitude oscillations. Numerically, only certain discretizations of the sphere can provide 
perfect quadrature formulae to compute the Fourier coefficients of band-limited functions on the 
sphere, sometimes with very efficient algorithms \274\ 275J. 

5.1.3 Spherical Scale-Space 

Similarly to what happened for signals or images, the first notion of "scale" on the sphere was 
imported from the Heat Dynamic that is also known on this space. In that framework, if a spherical 
function / E JL 2 (S 2 ) is considered the initial heat configuration, the spherical heat dynamics smooth 
it with time r > 0, conferring a scaling notion on this parameter. 

Interestingly, as for Euclidean spaces, the solution at time r > of the heat equation initialized 
to some function / E h 2 (S 2 ) is simply /(£,t) = ^ m /^ m (r)^ m (0, with ft m (r) = fime~ £ ^ 1)r 
and /(£, 0) = /(O- Alternatively, since for an axisymmetric filter h we have the spherical convolu- 
tion theorem 

the solution of the Heat Equation can also be obtained by a convolution by a specific kernel G°(0? 
coined spherical Gaussian of width y/r. It is defined in frequency by (G°)^ m = \J {21 + 1) /47T e -^+ 1 ) r . 

The link between the heat dynamics and the spherical convolution with the axisymmetric fil- 
ter G° T has been exploited by Biilow |276| to develop several specific spherical filters for feature 
detection, such as the Laplacian of Gaussian or the directional derivative of Gaussian. 
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5.1.4 Spectral Wavelets 

Freeden et al. \277\ I278j have fully exploited the connection between convolution and frequency 
filtering on the sphere to develop a continuous wavelet transform on the sphere. This is done 
by introducing a family of axisymmetric functions ^ a (€,), coined spherical wavelet, continuously 
indexed by a > 0, and such that J R+ \{^a)^\ 2 da/a = 1, (^a)oo = 0, Pl us additional regularity 
conditions. The wavelet coefficients of a function / E L, 2 (S 2 ) are then defined as Wf(a,£) = 
(f * il>a)(£)' The reconstruction is possible (almost everywhere) by 



In [27 7\ 1278] , an MRA on the sphere is also built by defining Quadrature Mirror Filters in the 
frequency domain. A spatial sub-sampling of the different subspaces of the MRA can also decrease 
the redundancy of the basis hence created. 

Following a similar approach, (isotropic) needlet frames introduced in \279\ \280\ 1281] represent 
another example of spectral wavelets, i.e., wavelets shaped in the Fourier domain. Needlets addi- 
tionaly offer relationships with quadrature formulae used to turn integrals of bandlimited functions 
into discrete summations. 

5.1.5 Stereographic Wavelets 

In the previous sections, the notion of scale in the processing of spherical data was always 
defined in the frequency domain, i.e., by dilating the frequency domain by a parameter, preventing 
a fine control of the spatial support of the filter. 

An alternative approach introduced by Antoine and Vandergheynst [282^ [283J defines the dila- 
tion directly in the spatial domain. The compactness of S 2 is respected, by introducing a stereo- 
graphic dilation. As illustrated on Fig. [30[(a) for point dilation, the stereographic dilation D a of a 
function g E h 2 (S 2 ) amounts to projecting g on the plane tangent at the North Pole by the stereo- 
graphic projection II, to applying there a Euclidean dilation d a by a scale a > 0, and to lifting the 
resulting function back to the sphere by II -1 [284]. Mathematically, [D a g](9, cp) = A(a, 9) g(0i/ a , ip), 
with tan 9 a /2 = a tan 9/2 and where A is a normalizing function such that ||D a ^||2 = ||p||. 

Given a mother wavelet ip E L 2 (S' 2 ) centered on the North pole, the proposed approach considers 
the joint action of translations, i.e., rotation operators R p in SO(3), and of the dilations D a on ip. 
The wavelet transform of / is therefore: 



where v is the Lebesgue measure on SO (3) and is a multiplicative operator function of ij) only 
and expressed in the Fourier domain [282]. For axisymmetric wavelets, this result simplifies by the 
fact that the action of R p on ip is controlled by two angles only. 




with (/) = ±f s2 f(£) d/z(£). 




l+cos<9 




z z 




(d) (e) 

Figure 30: (a) Stereographic dilation on S 2 . On the right, the (steerable) second directional 
derivative of Gaussian. The three images (b)-(d) are the basis elements, while the fourth in (e) is 
a linear combination of the first three yielding a rotation of 7r/4 around the North pole. 



Many wavelets may be defined on the sphere since it has been proved in [284] that any admissible 
wavelet on the plane L 2 (R 2 ) can be imported by inverse stereographic projection II -1 . A Laplacian 
of Gaussian (LoG), difference of Gaussians (DoG), Morlet Wavelet, and many other are generally 
used [2821 128511286) . Numerically, this spherical CWT is obtained thanks to the convolution theorem 
mentionned previously. This transform has been for instance intensively used in the analysis of the 
Cosmic Microwave Background (CMB), an astronomical signal remnant of some specific evolution 
phase of the Big Bang [283 128811289] . 

Wavelet frames can be developed in this theory by discretizing the scaling parameter a [285J. 
These frames, that do not subsample the spherical positions, have successfully served for the con- 
struction of invertible filter banks on the 2-Sphere [290] even if the stereographic dilation is not 
really compatible with the frequency description of the wavelets. 



5.1.6 Haar Transform on the Sphere 

The constructions of spherical wavelets described in the previous section make use of the Fourier 
decomposition on the sphere. It is possible to define wavelets directly over the spherical domain 



without Fourier analysis, using for instance the lifting scheme method [291J, see Sec. 5.3 This allows 
one to define spherical wavelets with a compact support, although the stability of the resulting 
transform is more difficult to control than over the planar domain. 

Inspired by this lifting scheme |291| , one can easily define a Haar basis on the sphere by consid- 
ering a family {M-^}j^j^o of embedded spherical triangulations that approximate a sphere [292J. 
These triangulations are obtained by a regular 1:4 refinement rule starting from an initial regular 
polyhedron A4°, and the edges are projected on the sphere to define spherical triangles. 

The corresponding spherical multiresolution defines Vj C L 2 (5 2 ) as the set of functions that 



are constant on each triangle of M 3 . Figure |3l| shows the linear projection of a spherical function 
on some of these multiresolution spaces. 
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j = 5 j = 4 j = 3 j = 2 

Figure 31: Projection on the spherical Haar multiresolution. 



Following the usual definition (Sec. 2.3.1), a Haar wavelet basis {ipj, n }n is an orthogonal basis 
of the detail space Wj such that Vj+i = Vj © Wj. The wavelet coefficients (^- n , /) are computed 
using a pyramid algorithm that mimics the usual Haar transform, except that for each triangle, 



one gathers three detail coefficients and one coarse scale coefficient. Figure [32] shows these Haar 
coefficients together with a comparison between spherical and planar non-linear approximations 
H T {f). 




Hr(f) (spherical) 



H T (f) (planar) 



Figure 32: Comparison of spherical and planar Haar approximations. The threshold T is adjusted 
so that Hx(f) is an approximation with a number of coefficients equal to 5% of the number of 
pixels in the high resolution planar image. 



5.1.7 Steerable Wavelets on the Sphere 

Finally, the sphere is compatible with the definition of steerable filters similarly to those defined 
in Sec. |3.2.2| for the plane. In particular, using the stereographic projection n introduced in the 
previous section, steerability on the sphere is also imported from the plane. This fact has been 
used in [284^ 293 , 294j to define differential and steerable filters useful to detect directional features 



in the Cosmic Microwave Background. An example of a steerable wavelet is given in Fig. 30'b-e). 
Spherical steerability may also be directly studied in the frequency domain with spectral dilation 
[295]. 
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5.1.8 Other Constructions 

It is impossible to cite the vast literature on multiscale decomposition on the sphere. Let us 
just quote some of them. Wavelets, ridgelets and curvelets have been translated on the sphere by 
Starck et al. |296j by using a particular spherical sampling, called HEALPix, locally similar to a 
square discretization. Locally supported biorthogonal wavelet bases have been also realized thanks 
to some radial projections of the planar faces of a cube on S 2 in [297J. 



5.2 Wavelets on General 2- Manifolds 

Given a two-dimensional manifold M, i.e., locally isomorphic to R 2 , authors in [298J describe 
how to define a Continuous Wavelet Transform (CWT) for function / : M —> C. 

Similarly to the way the stereographic dilation is defined for the sphere, the local dilation of 
a function around the point £ E M relies on the knowledge of a local and invert ible projection 

between A4 and its tangent plane T^M. on The desired dilation of scale a > therefore 
factorizes as £)(£, a ) — n^" 1 ^]!^ with d a the common Euclidean dilation of function in T^M. ~ R 2 . 

Given the Hilbert space H = L 2 (A4,d/x) of square integrable function on M, for a proper 
measure d/x, the CWT of a function / on M is then formally defined by correlating / with a set of 
prototype wavelets if;^ E % localized around any £ E M, i.e., 

Wf(M = (^a)J)n = [ d M (0/(0%a)(0> %a)=%,a)V>(0- 

JM 

The theoretical invertibility of this transform has however to be studied specifically in each case, 
i.e., given A4 and II^. Results exist for instance for the two-sheeted hyperboloid and the paraboloid 
in R 3 [299] . 



5.3 Lifting Scheme Wavelets on Meshed Surfaces 



The lifting scheme of Sweldens [300], described in Sec. 4.3, can be used to define wavelets on non- 
translation invariant geometries, including surfaces with complicated topologies. Lifted wavelets 
on surfaces are usually built on a semi-regular mesh grid, was first considered by Lounsbery et al. 
[301J, and then refined within the lifting framework by Schroder and Sweldens [291] . 

Semi-regular meshes {M J }j^j^o are obtained by a regular 1:4 refinement rule starting from an 
arbitrary control mesh A4°. Each edge of M.i is split into two sub-edges by vertex insertion to 
obtain the refined mesh ./VP -1 . The fine mesh A4 J is the sampling grid that stores the position of 



the surface points in space, and a signal / sampled at each grid point. Fig. [33j top row, shows an 
example of such a multiresolution mesh, obtained by a semi-regular remeshing of a high resolution 
input mesh. 



The lifting scheme described in Sec. 4.3 can be applied by storing the scaling coefficients aj on 
the grid point of the mesh Ai J , while the detail coefficients are stored on the complementary detail 
grid P J where A4 J_1 = M.^ UP J . The splitting of aj-i into a° and d° corresponds to assigning the 
values stored in .M J_1 to either M.i or VK The predict operator Pj used to compute the wavelet 
coefficients dj stored in P J is a local polynomial interpolator on a triangulation grid. The update 
operator Uj is computed by solving a linear system, to impose that moments of low orders, such as 
the mean, are preserved when moving from aj-i to aj. 
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This lifting wavelet transform computes the coefficients dj[n] = (^(j,™)? /) f° r a U scales < 
j < J and grid points n E T>K It corresponds to the projection of the signal / defined on 
the triangulated surface M.i onto a discrete biorthogonal wavelet frame B — {^(j,n)}j,n- These 
coefficients can be thresholded, and inverting the lifting steps creates an approximated signal /m 
with M non-zero coefficients. Although this approach works well in practice, the frame bounds of 
the resulting wavelet frame B are difficult to control, and Jm might be far from the best M-terms 
approximation. It is also difficult to guarantee the convergence of the wavelet atoms ^(j,™) to 
smooth functions, when J tends to — oo, and the mesh A4 J approximates a smooth surface. 

To perform surface approximation, one defines the signal aj at the finest scale as the position 
of the nodes on the surface. Each coefficient aj[n] E R 3 is thus a point in 3D space. The lifting 
transform can be applied to this vector- valued signal. Thresholding the resulting wavelet coefficients 



allows one to approximate the surface using few coefficients, as shown on Fig. [33j bottom row. If 
the lifting operators Pj and Uj do not depend on the position of the points on the surface, the 
resulting lifting wavelets can be used to perform 3D mesh compression |301} 1 291 j . 




(a) 



(b) 



(c) 



(d) 




(e) 



(f) 



(g) 



(h) 



Figure 33: Top row: example of semi-regular mesh {M J }j. Bottom row: example of surface 
approximation /m obtained by thresholding the lifted wavelets coefficients, where ./V is the number 
of vertices in M J . (a) M j ,j = -4. (b) M j ,j = -5. (c) MP, j = -6. (d) M j ,j = -7. (e) 
M/N = 100%. (f) M/N = 10%. (g) M/N = 5% (h) M/N = 2%. 
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5.4 Wavelets on Graphs 

Let us finally mention that wavelet transform has been extended to functions defined on the ver- 
tices of an arbitrary finite weighted graph. The latter may for instance generalize standard picture 
definition by describing two-dimensional pixel adjacencies. Maggioni et al introduced "diffusion 
wavelets" [302] , a general theory for wavelet decompositions based on compressed representations 
of powers of a diffusion operator such as the graph Laplacian. The constructed wavelet basis is 
made orthogonal by combining graph subsampling and Gram-Schmidt orthogonalization on each 
subsampled space. 

More recently, Hammond et al. [303] developed a general wavelet frame theory on such graphs 
thanks to the graph analogue of the Fourier domain, namely the spectral decomposition of the 
discrete graph Laplacian. Wavelets are defined in this frequency domain by dilating an "admissible" 
generating kernel. The final representation is redundant but wavelets can be shaped by changing the 
generating kernel. Moreover, for sparse graph Laplacian matrix, a fast wavelet transform avoiding 
the Laplacian spectral decomposition is developed. 

6 Conclusion 

A century after the discovery by Alfred Haar, and twenty years after the emergence of wavelets 
as genuine processing tools, major advances have been made in the improvement of natural images 
representations, aiming at enhanced understanding. 

Their common characteristic resides in uncovering multiscale and oriented features of natural 
images, through projections on a specific set of elongated atoms. The resulting dictionaries are 
thus often redundant, and may be coupled with sparsity enforcing priors, or adaptivity. They 
reveal a striking similarity with low level vision, where similar strategies are used to build powerful 
processing architectures. 

The availability of such a large number of transformations, that potentially extend the standard 
wavelet framework, leaves open the question of the best representation to process a given image. 
This choice is unfortunately data dependent, since the geometry of edges and textures varies sig- 
nificantly from natural to seismic or medical images. Selection of a representation, as well as its 
parameterization (number of scales, span of orientations, support in space or frequency), is also 
application dependent, and applications to inverse problems or pattern recognition typically im- 
pose strong design requirements on the dictionary. Their exhaustive comparison thus remains out 
of reach, with traditional methods from image processing or approximation theory only providing 
a partial answer. 

As a humble contribution to a subjective comparison, additional materials, full scale decompo- 
sition images, related links and associated toolboxes necessary to reproduce illustrations provided 
in this paper are available at [26]. Oddly enough, a common etymology of Szeged resides in an old 
Hungarian word for corner (szeg). At a turn in a wavelet century, A. Haar and F. Riesz might not 
have foreseen the harvest from their mathematical seeds. Image understanding is at the beginning 
of reaping their fruits. 
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